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Using high-resolution data from the Galactic Arecibo L-Band Feed Array HI (GALFA-Hl) survey, 
we show that linear structure in Galactic neutral hydrogen (Hi) correlates with the magnetic field 
orientation implied by Planck 353 GHz polarized dust emission. The structure of the neutral 
interstellar medium is more tightly coupled to the magnetic field than previously known. At high 
Galactic latitudes, where the Planck data are noise-dominated, the Hi data provide an independent 
constraint on the Galactic magnetic field orientation, and hence the local dust polarization angle. We 
detect strong cross-correlations between template maps constructed from estimates of dust intensity 
combined with either Hi-derived angles, starlight polarization angles, or Planck 353 GHz angles. 

The Hi data thus provide a new tool in the search for inflationary gravitational wave H-mode 
polarization in the cosmic microwave background, which is currently limited by dust foreground 
contamination. 


The cosmic microwave background (CMB) is the per¬ 
vasive residual radiation from the formation of the Uni¬ 
verse. The detection of primordial 5-mode polarization 
in the CMB is a major goal of contemporary cosmology. 
This signal is imprinted at the surface of last scatter¬ 
ing by perturbations from gravitational waves generated 
during the epoch of inflation, a period of rapid expansion 
in the early Universe [1-3]. An inflationary gravitational 
wave (IGW) 5-mode measurement would be the first di¬ 
rect evidence of inflation. A number of experiments are 
pursuing the signal, using ground-based (e.g. ABS [4]; 
Advanced ACT [5]; BICEF2/Keck Array [6]; CLASS [7]; 
POLARBEAR [8]; SPT-3G [9]), balloon-borne (EBEX 
[10]; SPIDER [11]), and space telescopes {Planck [12]). 

Unfortunately, our view of the polarized CMB is ob¬ 
scured by contaminating foregrounds. Eor IGW 5-mode 
searches at frequencies ^ 100 GHz, the largest fore¬ 
ground is Galactic polarized dust emission. Aspherical 
dust grains in the Milky Way align their short axes with 
the ambient magnetic field, and interstellar radiation is 
absorbed and reradiated by the dust as partially polar¬ 
ized light. The BICEP2 collaboration claimed a measure¬ 
ment of primordial 5-modes [13], but subsequent analy¬ 
ses determined that the detection could be attributed en¬ 
tirely to Galactic dust [14, 15]. A detailed understanding 
of the foreground polarization signal is required before a 
definitive IGW 5-mode detection can be achieved. Pur¬ 
suant to that goal, the Planck satellite recently mapped 
the full sky at 353 GHz, a frequency dominated by ther¬ 
mal dust emission. These data can be used to subtract 
the foreground polarization pattern from lower-frequency 
CMB observations. To optimize the chance of primordial 
5-mode detection, experiments should target the “clean¬ 
est” regions of sky: areas where there is relatively little 
polarized dust, and where the dust polarization structure 
is measured with high signal-to-noise. The Planck maps 
are limited in this regard, because the Planck polarized 


signal is noise-limited at high Galactic latitudes, where 
the dust column is lowest. Thus IGW 5-mode searches 
are plagued by a trade-off: the regions of lowest fore¬ 
ground amplitude are also the regions with the poorest 
foreground constraints. 

We present an entirely new method for constraining 
Galactic foregrounds. Using only the morphology of dif¬ 
fuse neutral hydrogen (Hi) structures, we predict the ori¬ 
entation of polarized dust emission at high precision over 
a range of angular scales. In parallel with existing mea¬ 
surements of polarized CMB foregrounds, our recovery of 
the dust polarization angle will increase the precision of 
foreground models. This is especially valuable in regions 
where the Planck 353 GHz data are noise-limited. 

This work follows the discovery that linear structures 
in Hi are elongated in the direction of the interstellar 
magnetic field as probed by starlight polarization [16]. 
Here, we demonstrate that Hi orientation is well corre¬ 
lated with the Planck 353 GHz polarization angle across 
a region of high Galactic latitude sky. Note that the 
Planck data enable quantitative conclusions beyond the 
previous work, which considered only 153 sparsely sam¬ 
pled starlight polarization measures over 1,300 deg^ of 
sky. Also, polarized dust emission samples the full line 
of sight, whereas starlight polarization only traces the 
magnetic field out to the distance of the star. The rela¬ 
tionship between dust and Hi in the interstellar medium 
(ISM) is deeper than their correlation in column density 
(e.g. [17]), which is already used to estimate the ampli¬ 
tude of polarized dust emission [14, 18]. Small dust grains 
and long plumes of Hi are both aligned by the magnetic 
field, though the mechanism for aligned Hi structure for¬ 
mation is not yet well understood. 

The slender linear features that best trace the orien¬ 
tation of the Galactic magnetic field are only revealed 
by high spatial and spectral resolution Hi maps (see [16] 
for details). We use data from the Galactic Arecibo L- 
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FIG. 1. Plane-of-sky magnetic field orientation as predicted by Orht (top) and ^353 (bottom). Color maps are integrated HI 
column density from v — —61.5 km s“^ to +61.5 km s“^ {Nhi [cm“^]), and dust opacity (t353)- Planck and RHT Q and U 
values are smoothed with a FWHM = 1° Gaussian kernel, then used to construct ^353 and Orht, which are visualized using 
line integral convolution (LIC; [20]). The high latitude {b > 70) behavior of the ^353 LIC pattern is due to Planck noise. White 
pseudovectors represent starlight polarization angles. Galactic latitude lines lie at 5 = 30°, 50°, 70°, from left to right. Galactic 
longitude lines lie at I = 80°, 50°, 20°, from top to bottom. 

Band Feed Array Hi survey (GALFA-Hl; [19]) with the 
Arecibo 305m radio antenna. GALFA-Hl has an angular 
resolution of FWHM 4', a spectral resolution of 0.18 
km s“^, and a brightness temperature noise of ^ 140 mK 
rms per 1 km s“^ integrated channel over 13,000 deg^ of 
sky. This work uses data from the forthcoming second 
data release [ 21 ]. 

We analyze 353 GHz polarization data obtained by the 
Planck satellite’s High Frequency Instrument (HFI) [18]. 

These data have an angular resolution of FWHM 5', 
comparable to GALFA-Hl. We transform the Planck 
data from Galactic to Equatorial coordinates [ 22 ]. For 
all analyses, we apply a mask constructed from the union 
of all point source masks provided for each HFI channel 
in both temperature and polarization. 

We additionally consider 126 optical starlight polariza¬ 
tion measures in this region [23]. Starlight is polarized 
parallel to the magnetic field by the preferential absorp¬ 
tion of aligned grains. 

We quantify the orientation of GALFA-Hl structures 
using the Rolling Hough Transform (RHT), a machine 
vision technique [16]. The RHT runs on image data, and 
outputs R{0), linear intensity as a function of angle, for 
every pixel in the input map. For a detailed description 
of the RHT we refer the reader to [16]. 

For this work we select a 1,278 deg^ region of the 


GALFA-Hl sky. The region, which spans right ascension 
195° to 265° and declination 19.1° to 38.3°, stretches 
from b = 30° above the Galactic plane to b = 81.7°, 
nearly Galactic zenith. We analyze this GALFA-Hl re¬ 
gion from —13.5 km s“^ to +13.5 km s“^, binned in 3.0 
km s“^ integrated velocity channels. 

Linear polarization data can be fully described by ei¬ 
ther a polarization angle -0 and polarized intensity P 
or by the Stokes parameters Q and 77, where -0 = 
1/2 arctan(7//(5) and P U^. We define from 

the RHT output 


Qrht = ^ 

1 cos (26») • R (0) dd 


Urht = 

f sin (20) ■ R (0) d0, 

( 1 ) 


where values are calculated for each point in the image 
data. We process each velocity channel with the RHT 
and add the resulting Qrht and Urht maps. 

We define Orht = ^^^^^^{Urht/Qrht)^ an es¬ 

timate for the orientation of the magnetic field derived 
solely from Hi data. We compare this value to 6 ^ 353 , a 90° 
rotation of the polarization angle obtained from Q 353 and 
^353 (we use the lAU polarization definition). The polar¬ 
ization angle of dust emission is conventionally taken to 
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FIG. 2. Normalized histograms of SO = ^353 — Orht in 1° 
bins at resolutions of FWHM = 5', lb', and 30b The Gaus¬ 
sian ht to the FWHM =15^ histogram shown has a standard 
deviation a = 19.4°. 

be 90° from the orientation of the local Galactic magnetic 
field (however, see [24] and references therein). 

We calculate ^353 and Orht for the region described. 
Figure 1 shows a map of each of these quantities on the 
sky, along with starlight polarization angles. Although 
derived from independent data, these three estimates for 
the plane-of-sky magnetic field orientation trace one an¬ 
other remarkably well. Figure 2 shows histograms of 
SO = ^353 — Orht- We construct ^353 and Orht from Q 
and U maps smoothed to three different resolutions. For 
Gaussian smoothing kernels of FWHM = 5', 15', and 30', 
we find Gaussian fits to the SO histogram with standard 
deviation a = 30.2°, a = 19.4°, and a = 14.4°, respec¬ 
tively. We run a Monte Carlo analysis to determine the 
pure Planck noise contribution to SO, and find this noise 
is responsible for a Gaussian component with a = 16.0°, 
a = 6.1°, and a = 3.5° for each respective smoothing 
kernel. Thus as the data are smoothed to larger angular 
scales, ^353 and Orht obtain ever better agreement, and 
a non-negligible fraction of the SO scatter is solely due 
to Q 353 and 1/353 measurement noise. The SO histograms 
are centered at about —3° to —4°. This small offset from 
zero may be due to either residual systematics in the 353 
GHz map [25] or true systematic differences between ^353 
and Orht- 

To further characterize the relationship between RHT, 
Planck, and starlight polarization angles, we construct 
simple template maps and compute cross-power spectra 
between them. We construct the templates using the 
Planck 353 GHz intensity, 1353 . A full polarization tem¬ 
plate would also require an estimate of the polarization 
fraction, p = P/1, but since our goal is to isolate the 
polarization angle information, we set p = 1 in all tem¬ 
plates. (Over a small patch of sky, p ^ constant is a 
reasonable approximation, and one can simply re-scale 
our power spectra for a given value of p.) Furthermore, 
measuring P from the Planck data is non-trivial, as sim¬ 


ple estimators are noise-biased (e.g. [26]). The templates 
are 

Q = ^353 cos(2'0) and U = 1353 sin(2'0), ( 2 ) 

where 'ip is either the RHT, Planck, or starlight polar¬ 
ization angle. For all templates, we smooth the Q and 
U data to a common resolution of FWHM = 4° be¬ 
fore computing p;. This prohibits small-scale noise in 
Q and U from contaminating the templates on large 
scales via the harmonic-space convolution implied by the 
real-space map multiplication in Eq. (2). To avoid noise 
biases, we measure cross-correlations between templates 
constructed from independent half-mission splits of the 
Planck data. 

We apply a common mask to all template maps, con¬ 
sisting of the Planck point source mask and a mask re¬ 
moving regions that are more than 7° from starlight data, 
regions where the integrated RHT intensity is zero, and 
the edges of the region, where RHT artifacts could arise. 
The total unmasked sky subtends 1,181 deg^, 92% of the 
original area. We apodize the mask with a Gaussian 
taper of FWHM = 15'. We use polspice [27] to com¬ 
pute EE and BB power spectra {Ci, where multipole 
i is the harmonic variable conjugate to angular scale), 
corresponding to the usual curl-free and divergence-free 
decompositions of polarization data [1, 3], respectively. 
We calibrate the polspice internal parameters using 100 
simulations of polarized dust power spectra with proper¬ 
ties matching recent Planck measurements [18, 28]. We 
bin the measured power spectra in four logarithmically 
spaced multipole bins between ^ = 40 and i = 600 (cen¬ 
tered at = 59,116,229, and 451). Error bars are cal¬ 
culated in the Gaussian approximation from the auto¬ 
power spectra of the template maps used in each cross¬ 
correlation. Sample variance is not included in the error 
bars, as our interest is in comparing measurements of the 
same modes on the sky. 

Eigure 3 shows cross-power spectra for the template 
maps constructed from RHT, Planck, and starlight data. 
We refrain from fitting a model to the data, as we 
have not considered p in our templates, but instead con¬ 
sider the relative amplitudes of the cross-power spectra. 
Eor the Planck-only templates and the Planck-KWT, 
Planck-stailight, and RHT-starlight cross-correlations, 
respectively, we detect the E-mode power spectrum at 
70cr, 55cr, 40(7, and 40(7 significance. We detect the B- 
mode power spectrum at 65(7, 60(7, 50(7, and 40(7 sig¬ 
nificance. We verify that template maps constructed 
with random angles yield a cross-power spectrum con¬ 
sistent with zero (even when using the true 1353 data in 
the random-angle templates). We compare the template 
cross-power spectra with the actual EE and BB power 
spectra measured directly from Qsss and t /353 and infer 
a mean p ^ 5%, which is reasonable for this region [18]. 

Although significant cross-correlations are detected for 
all templates in Eigure 3, the Planck-only templates 
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FIG. 3. Cross-power spectra of polarization template maps 
constructed from 7353 and either Planck (' 0353 ), RHT {'iPrht), 
or starlight polarization data (Eq. 2). Shown are E-mode 
(circles) and E-mode (squares) components. Significant (40- 
70(7) cross-correlations are detected in all cases. 


yield higher amplitudes than the cross-correlations with 
RHT- or starlight-based templates. While this could be 
due to physical differences between angles, we note that 
the RWT-Planck and RHT-starlight cross-correlations 
yield similar results (especially at low-7), suggesting that 
the Planck-only templates’ power spectra could be sys¬ 
tematically biased. Because the angle construction re¬ 
lies on the 7 / 353 /Q 353 ratio, it is sensitive to any ef¬ 
fect that modifies the zero point of the maps. Such ef¬ 
fects could include gain calibration drifts or intensity- 
to-polarization leakage that varies over the sky, both of 
which are known to be present in the Planck data at 
some level [25]. Indeed, scan-synchronous systematics 
have been detected in Planck temperature data [29, 30], 
and maps of SO = ^353 — Orht present clear visual ev¬ 
idence of residuals that are highly correlated with the 
Planck scan directions. We leave a detailed consider¬ 
ation of these systematic effects on the Planck angles 
to future work. Note that direct measurements of EE 
and BB power spectra from Q 353 and 7/353 are more im¬ 
mune to these systematics than the angle construction, 
but we require the latter method to compare Planck data 
in a straightforward way to the RHT- and starlight-based 
templates. 

The KRT-Planck cross-power spectra yield an ampli¬ 
tude ratio 2 , a result consistent with the 

Planck 353 GHz measurement [28] , though this must be 
interpreted with caution as we have not modeled p in 
our templates. Many current models of the dust polar¬ 
ized sky [31, 32] predict equal E- and E-mode amplitudes 
[28]. Hi orientation preserves the nonunity EE/BB ra¬ 
tio, suggesting that ISM structure is a crucial missing 
component of these models. The preferential alignment 
of Planck filamentary dust structures with the magnetic 
field [33] supports this conclusion. Our work underscores 



FIG. 4. Angle uncertainties averaged over 2 ° Galactic latitude 
bins, normalized by their respective median values (dashed 
line). 


the need for a deeper understanding of the interplay be¬ 
tween ISM phenomena and polarized dust. 

IGW E-mode experiments often target the high Galac¬ 
tic latitude sky, where Planck data cannot distinguish be¬ 
tween the most promising potential targets [34]. Figure 4 
shows the relative Galactic latitude dependence of uncer¬ 
tainties in 6^353 and Orht^ where the Orht uncertainty 
is propagated from the variance in R (0). With sensitive 
measurements at high latitudes, Orht maps can be used 
to assess the structure of the magnetic field in targeted 
regions of sky. 

Our results indicate that full foreground templates 
with higher signal-to-noise than the (^353 and 7/353 maps 
can be constructed by combining Orht with other data 
describing P. A scale-dependent modeling of p and I 
from a combination of 1353 , P 353 , and Hi data may en¬ 
able the extension of this work to full polarized dust fore¬ 
ground maps. Such templates should remove CMB and 
cosmic infrared background emission from 7353 , which we 
neglect here. We can also replace 7353 in Eq. 2 with an 
unbiased estimator of P 353 (e.g. [26, 35]). P is theoret¬ 
ically determined by the dust column along the line of 
sight, traced by 7, and the tangledness of the magnetic 
field along the line of sight, where more tangled fields 
cause greater depolarization. Nhi is a powerful proxy for 
^ 353 , particularly at high Galactic latitudes where dust 
emission is low and the expected depletion of Hi into a 
molecular state is minimal. Changes in Orht for differ¬ 
ent Hi velocity channels may indicate line-of-sight field 
tangling, and may elucidate the physical origin of varia¬ 
tions in p by isolating components of the magnetic field. 
This will be the subject of future work, and may lead to 
further Hi constraints on P353- 

In this work we demonstrate that Hi orientation cor¬ 
relates with Planck 353 GHz polarization angle. We will 
process Hi data from the full Arecibo sky in a forth¬ 
coming work, as it overlaps with several CMB experi¬ 
ments. Lower resolution Hi surveys such as GASS [36] 
and EBHIS [37] can be used on other regions of the sky, 
although they do not trace the Galactic magnetic field 
as precisely as the high resolution GALFA-Hl data [16]. 
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Soon, Galactic all-sky maps from Square Kilometer Ar¬ 
ray pathfinders [38] will be ideal for Hl-based foreground 
maps. 
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SUPPLEMENTAL MATERIAL 
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FIG. 5. Approximate observing regions of various CMB B-mode experiments in Equatorial coordinates. Background image is 
the log of the Planck 353 GHz dust intensity [K]. The GALFA-Hl full-sky region is overlaid. The Effelsberg-Bonn Hi Survey 
(EBHIS) and the Galactic All-Sky Survey (GASS) cover the entire Northern and Southern Equatorial skies, respectively. 


In this work we run the RHT using an unsharp mask kernel diameter Dx = 15', a rolling window size Dw = 75', 
and an intensity threshold Z = 70% (see [16]). Under variation of both the velocity channel binning as well as the 
RHT parameters, our results remain qualitatively unchanged. 

Uncertainties for P and p are defined by 


crp = —^ 


( 3 ) 


and 




28.65° 


Q‘^cfuu‘^ + U‘^(Jqq‘^ 
+ U‘^auu^ 


crp 


( 4 ) 


where we neglect QU covariance. We compute ctqq and from half-mission splits of the Planck data, following 
the procedure outlined in [25]. The analogous quantities cfqq^ and are computed from the variance in the 

RHT spectrum, as 


= jcos\2e) ■ R{e) de (5) 

and 

= j sin^ (26») • R (6») dO. (6) 










































We define the difference between 6^353 and Orht as 


SO = -arctan 
2 


sin( 26 > 353 ) cos{ 20 rht) - cos( 26 > 353 ) sm{ 20 RHT) 

COs(20353) C0s{2diiHT) + sin(20353) sin{2diiHT) 


This equation properly accounts for the 180° degeneracy in polarization angle. 



